% 绘制函数pip的图像
x = linspace(0,3,100);
y = pip(x);
figure(1);
plot(x, y)
xlabel('x');
ylabel('pip(x)');
text(0.5, pip(0.5), num2str(pip(0.5)), 'color','g');
text(1, pip(1), num2str(pip(1)), 'color','g');
title('土壤温度随深度的变化图');
saveas(1, '土壤温度随深度的变化图.jpg');

% 计算温度为0时埋管深度
a = fzero(@pip, [0.5, 1]);
fprintf('埋管深度为%.2f米时水管刚好不冻结\n', a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数功能：  插值求导
% 参数意义：  x: 水管埋藏的深度
% 返回值：    y：水管在深度h对应的温度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = pip(x)
% 热传导系数
alpha = 0.138e-6;
% 正常土壤温度
T_i = 20;
% 寒冷季节地面温度
T_s = -15;
% 寒流持续时间（单位为秒）
t = 60 * 24 * 60 * 60;
% 深度h对应的土壤温度计算公式
y = T_s + (T_i - T_s) * erf(x / (2 * sqrt(alpha * t)));
end

